      MODULE INOUT_PARAMETERS
      
      integer(4), parameter :: lake_subr_unit_min = 11001
      integer(4), parameter :: lake_subr_unit_max = 12000
      integer(4), parameter :: lake_mon_out_unit_min = 12001
      integer(4), parameter :: lake_mon_out_unit_max = 13000
      integer(4), parameter :: lake_day_out_unit_min = 13001
      integer(4), parameter :: lake_day_out_unit_max = 14000
      integer(4), parameter :: lake_hour_out_unit_min = 14001
      integer(4), parameter :: lake_hour_out_unit_max = 15000
      integer(4), parameter :: lake_everystep_out_unit_min = 15001
      integer(4), parameter :: lake_everystep_out_unit_max = 16000
      integer(4), parameter :: lake_series_out_unit_min = 16001
      integer(4), parameter :: lake_series_out_unit_max = 17000
      
      END MODULE INOUT_PARAMETERS


      MODULE PHYS_PARAMETERS
      
      real(8) :: whc, hcond, sncr
      
      SAVE

      contains
      SUBROUTINE DEFINE_PHYS_PARAMETERS()
      implicit none

!cccccccccccccccccccccc SNOW cccccccccccccccccccccccc
      whc = 0.04               ! water holding capacity of snow
      hcond = 0.01             ! hydraulic conductivity in snow (m/s)
      SNCR = 0.01              ! M OF WATER, FROM WHICH SNOW COVERS ALL THE CELL

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      END SUBROUTINE DEFINE_PHYS_PARAMETERS
      END MODULE PHYS_PARAMETERS

 
      MODULE PHYS_CONSTANTS
      
      real(8) :: g
      real(8) :: omega
      real(8) :: roa0, row0, roi
      real(8) :: cw, ci, csnow
      real(8) :: Lwi, Lwv, Liv
      real(8) :: kappa
      real(8) :: roughnessofwater
      real(8) :: sigma
      real(8) :: cp
      real(8) :: aMag, bMag
      real(8) :: Rd, Rwv, R_univ
      real(8) :: lami, lamw0
      real(8) :: niu_atm
      real(8) :: roughness0
      real(8) :: albedoofice, albedoofsnow, albedoofwater, albedoofsoil
      real(8) :: emissivityofice, emissivityofwater 
      real(8) :: emissivityofsnow, emissivityofsoil
      real(8) :: aMagw, bMagw
      real(8) :: aMagi, bMagi
      real(8) :: alsal, almeth, aloxyg, alcarbdi
      real(8) :: sabs
      real(8) :: Kelvin0
      real(8) :: pref
      real(8) :: esatsurf0
      
!     Derived constants
      real(8) :: ci_m_roi, cw_m_row0
      real(8) :: row0_m_Lwi, roi_m_Lwi
      real(8) :: row0_d_roi, roi_d_row0
      real(8) :: roa0_d_row0, roa0_d_roi
      real(8) :: row0_d_roa0
      real(8) :: Rd_d_cp, Rd_d_Rwv
      
      
      SAVE
      
      contains
      SUBROUTINE DEFINE_PHYS_CONSTANTS 
      implicit none

      g     = 9.814    ! acceleration due to gravity,               m/s**2
      omega = 7.29d-5  ! angular velocity of the Earth's rotation,  1/s
      sigma = 5.67d-8  ! Stefan-Boltzman constant,                  W/(m**2*K**4)

      roa0  = 1.273    ! reference density of air,                  kg/m**3
      row0  = 1000.    ! reference density of water,                kg/m**3
      roi   = 917.     ! density of ice,                            kg/m**3

      cw    = 3990.    ! heat capacity of water,                    J/(kg*K)
      ci    = 2150.    ! heat capacity of ice,                      J/(kg*K)
      cp    = 1005.    ! heat capacity of air at constant pressure, J/(kg*K)
      csnow = 2100.    ! the snow specific heat content,            J/(kg K)

      lami  = 2.2      ! molecular conductivity of ice,             J/(m*s*K)
      lamw0 = 0.561    ! molecular conductivity of water,           J/(m*s*K)
      niu_atm = 1.8d-5 ! molecular viscosity of air,                m**2/s

      Lwi   = 267900.  ! latent heat of freezing,                   J/kg
      Lwv   = 2501000. ! latent heat of evaporation,                J/kg
      Liv   = 2834600. ! latent heat of sublimation,                J/kg

      kappa = 0.38     ! Karman constant,                           n/d 
      Rd    = 287.05   ! gas constant for dry air,                  J/(kg*K)
      Rwv   = 461.91   ! gas constant for water vapor,              J/(kg*K)
      R_univ = 8.31441 ! universal gas constant,                    J/(mol*K)

      aMagw = 7.6326   ! coefficient for Magnus formula for water surface
      bMagw = 241.9    ! coefficient for Magnus formula for water surface
      aMagi = 9.5      ! coefficient for Magnus formula for ice   surface
      bMagi = 265.5    ! coefficient for Magnus formula for ice   surface

      albedoofwater = 0.07 ! albedo of water surface,               n/d ! 0.07
      albedoofice   = 0.60 ! 0.5! albedo of ice   surface,               n/d ! 0.35
      albedoofsnow  = 0.60 ! 0.8! albedo of snow  surface,               n/d ! 0.62
      albedoofsoil  = 0.3  ! albedo of soil  surface,               n/d

      emissivityofwater = 0.99 ! emissivity of water surface,       n/d ! 0.95
      emissivityofice   = 0.99 ! emissivity of ice   surface,       n/d ! 0.95
      emissivityofsnow  = 0.99 ! emissivity of snow  surface,       n/d
      emissivityofsoil  = 0.9  ! emissivity of soil  surface,       n/d
      
      roughness0 = 1.d-3 ! the reference roughness of water, m
      
      Kelvin0 = 273.15
      pref    = 1.d+5 ! Reference atmospheric pressure, Pa
      esatsurf0 = 610.7 ! Saturation partial water vapor pressure at 0 Celsius, Pa

      sabs         = 0.1 !0.1 ! the part of solar radiation, absorbed at the water surface
      alsal        = 1.  ! the ratio of eddy diffusivity for salinity to those of heat
      almeth       = 1.  ! the ratio of eddy diffusivity for dissolved methane to those of heat
      aloxyg       = 1.  ! the ratio of eddy diffusivity for dissolved oxygen to those of heat
      alcarbdi     = 1.  ! the ratio of eddy diffusivity for dissolved carbon dioxode to those of heat
      
      ci_m_roi  = ci * roi
      cw_m_row0 = cw * row0
      row0_m_Lwi = row0 * Lwi
      roi_m_Lwi = roi * Lwi
      
      row0_d_roi = row0/roi
      roi_d_row0 = roi/row0
      roa0_d_row0 = roa0/row0
      roa0_d_roi = roa0/roi
      row0_d_roa0 = row0/roa0
      
      Rd_d_cp = Rd/cp
      Rd_d_Rwv = Rd/Rwv

      END SUBROUTINE DEFINE_PHYS_CONSTANTS
      END MODULE PHYS_CONSTANTS


      MODULE NUMERIC_PARAMS
       integer(4), PARAMETER:: KL = 1,ML = 131,MS = 100,NT = 52592, &
       & Num_Soil = 11,Num_Veget = 13 
       real(8), parameter:: dznorm = 0.05, dzmin = 0.01, &
       & UpperLayer = 0.01 ! UpperLayer = 0.1, UpperLayer = 0.008 
       
       real(8), parameter :: min_ice_thick = 0.01 
       real(8), parameter :: min_water_thick = 0.01
       real(8), parameter :: T_phase_threshold = 1.d-5

!         ML --- total number of layers
!         MS --- maximal number of layers in snow
!         dznorm --- standard layer depth in snow (m)
!         dzmin --- min layer depth in snow (m)
!         depth --- depth of the lowest layer in soil (m)
!         UpperLayer --- upper soil layer depth (m)

      END MODULE NUMERIC_PARAMS
      

      MODULE ATMOS

      real(8) :: tempair
      real(8) :: relhum
      real(8) :: humair
      real(8) :: pressure
      real(8) :: hw
      real(8) :: xlew
      real(8) :: cdmw
      real(8) :: uwind
      real(8) :: vwind
      real(8) :: wind, wind10
      real(8) :: longwave
      real(8) :: Radbal
      real(8) :: hflux
      real(8) :: surfrad
      real(8) :: hbal
      real(8) :: zref
      real(8) :: botflux
      real(8) :: cdmw2
      real(8) :: shortwave
      real(8) :: precip
      real(8) :: Sflux0
      real(8) :: Elatent
      real(8) :: Radbal_surf
      real(8) :: eflux
      real(8) :: eflux0
      real(8) :: eflux0_kinem
      real(8) :: hskin
      real(8) :: tau
      real(8) :: velfrict
      real(8) :: velfrict_prev
      real(8) :: cdm1
      real(8) :: turb_density_flux
     
      SAVE 

      END MODULE ATMOS


      MODULE DRIVING_PARAMS
      
      use INOUT_PARAMETERS, only : &
      & lake_subr_unit_min, &
      & lake_subr_unit_max
      
! Switches: 1. PBL parameterization
!              PBLpar =1 (Businger-Dayer formulas (Monin-Obukhov theory) for exchange coefficients)
!              PBLpar =11(Businger-Dayer formulas with shallow water correction after Panin(19..))      
!              PBLpar =2 (formulation from NH3d)
!              PBLpar =21(formulation from NH3d with shallow water correction after Panin(19..))
!           2. Relative to water currents wind
!              relwind =1 (relative wind is off)
!              relwind =2 (relative wind is on)
!           3. Turbulent mixing parameterization
!              Turbpar =1 (analytyical profile of turbulent conductivity)
!              Turbpar =2 ("E-eps" parameterization of turbulent conductivity)        
!              Turbpar =3 : Nickuradze (NICK) formulation: Rodi (1993)    
!              Turbpar =4 : Parabolic (PARAB) formulation: Engelund (1976)
!              Turbpar =7 : RNG (re-normalization group) formulation: Simoes (1998)     
    
      real(8) :: dt_out ; logical :: ok_dt_out = .false.
      real(8) :: sensflux0 ; logical :: ok_sensflux0 = .false.
      real(8) :: momflux0 ; logical :: ok_momflux0 = .false.
      real(8) :: kwe ; logical :: ok_kwe = .false.
      real(8) :: d_surf ; logical :: ok_d_surf = .false.
      real(8) :: d_bot ; logical :: ok_d_bot = .false.
      real(8) :: depth ; logical :: ok_depth = .false.

!     The group of tributaries characteristics

      integer(4) :: N_tribin !; logical :: ok_N_tribin = .false.
      integer(4) :: N_tribout !; logical :: ok_N_tribout = .false.
      integer(4) :: tribheat ; logical :: ok_tribheat = .false.
      real(8), allocatable :: &
      & U_tribin     (:,:), &
      & U_tribout    (:,:), &
      & T_tribin     (:,:), &
      & width_tribin (:,:), &
      & width_tribout(:,:)  
      real(8), allocatable :: zTinitprof(:)
      real(8), allocatable :: Tinitprof(:)
      real(8), allocatable :: zgrid_out(:)
      real(8), allocatable :: zgridsoil_out(:)

!     In perspective, it is expected, that area_lake be an array,
!     as the horizontal section of a lake varies with depth

      integer(4) :: nunit = lake_subr_unit_min
      
      integer(4) :: runmode ; logical :: ok_runmode = .false.
      integer(4) :: PBLpar ; logical :: ok_PBLpar = .false.
      integer(4) :: waveenh ; logical :: ok_waveenh = .false.
      integer(4) :: relwind ; logical :: ok_relwind = .false.
      integer(4) :: Turbpar ; logical :: ok_Turbpar = .false.
      integer(4) :: stabfunc ; logical :: ok_stabfunc = .false.
      integer(4) :: kepsbc ; logical :: ok_kepsbc = .false.
      integer(4) :: varalb ; logical :: ok_varalb = .false.
      integer(4) :: skin ; logical :: ok_skin = .false.
      integer(4) :: massflux ; logical :: ok_massflux = .false.
      integer(4) :: ifrad ; logical :: ok_ifrad = .false.
      integer(4) :: sedim ; logical :: ok_sedim = .false.
      integer(4) :: M ; logical :: ok_M = .false.
      integer(4) :: Mice ; logical :: ok_Mice = .false.
      integer(4) :: ns ; logical :: ok_ns = .false.
      integer(4) :: nstep_keps ; logical :: ok_nstep_keps = .false. ! The number of timesteps of k-epsilon parmeterization per on model timestep
      integer(4) :: nscreen ; logical :: ok_nscreen = .false.
      integer(4) :: SoilType ; logical :: ok_SoilType = .false.
      integer(4) :: Tinitlength  ; logical :: ok_Tinitlength = .false.
      integer(4) :: dyn_pgrad  ; logical :: ok_dyn_pgrad = .false.

      integer(4) :: monthly_out ; logical :: ok_monthly_out = .false.
      integer(4) :: daily_out ; logical :: ok_daily_out = .false.
      integer(4) :: hourly_out ; logical :: ok_hourly_out = .false.
      integer(4) :: time_series ; logical :: ok_time_series = .false.
      integer(4) :: everystep ; logical :: ok_everystep = .false.
      integer(4) :: turb_out ; logical :: ok_turb_out = .false.
      integer(4) :: scale_output ; logical :: ok_scale_output = .false.
      integer(4) :: ngrid_out ; logical :: ok_ngrid_out = .false.
      integer(4) :: ngridsoil_out ; logical :: ok_ngridsoil_out = .false.
      integer(4) :: assim ; logical :: ok_assim = .false.
      integer(4) :: as_window ; logical :: ok_as_window = .false.
      integer(4) :: error_cov ; logical :: ok_error_cov = .false.
      integer(4) :: zero_model ; logical :: ok_zero_model = .false.
      integer(4) :: year_accum_begin, year_accum_end
      integer(4) :: month_accum_begin, month_accum_end
      integer(4) :: day_accum_begin, day_accum_end
      integer(4) :: hour_accum_begin, hour_accum_end
      logical :: ok_accum_begin, ok_accum_end
      
      logical :: all_par_present = .true.

      character(len=40) :: path ; logical :: ok_path = .false.
      character(len=60) :: setupfile ! ; logical :: ok_setupfile = .false.
      
      real(8),    external :: getvarval
      integer(4), external :: igetvarval

      SAVE
      
      contains
      SUBROUTINE DEFINE_DRIVING_PARAMS()
      
!     DEFINE_driving_params reads files with driving parameters

      implicit none
      character(len=200) :: line
      logical :: firstcall
      data firstcall, line/.true., 'begin file'/
      
      if (firstcall) then
       call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,nunit)
       open  (nunit,file='setup_lake.dat',status='old')
       read  (nunit,'(a)') line
       read  (nunit,'(a)') setupfile
       close (nunit)
       open  (nunit,file=setupfile,status='old')
       do while (line /= 'end')
         read (nunit,'(a)') line
         call READPAR(line)
       enddo
       close (nunit)
       
       if (.not.ok_dt_out) then
         write(*,*) 'The parameter dt_out is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_sensflux0) then
         write(*,*) 'The parameter sensflux0 is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_momflux0) then
         write(*,*) 'The parameter momflux0 is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_kwe) then
         write(*,*) 'The parameter kwe is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_d_surf) then
         write(*,*) 'The parameter d_surf is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_d_bot) then
         write(*,*) 'The parameter d_bot is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_depth) then
         write(*,*) 'The parameter soil_depth is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_tribheat) then
         write(*,*) 'The parameter tribheat is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_runmode) then
         write(*,*) 'The parameter runmode is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_PBLpar) then
         write(*,*) 'The parameter PBLpar is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_waveenh) then
         write(*,*) 'The parameter waveenh is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_relwind) then
         write(*,*) 'The parameter relwind is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_Turbpar) then
         write(*,*) 'The parameter Turbpar is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_stabfunc) then
         write(*,*) 'The parameter stabfunc is missing in setup file'
         all_par_present = .false.
       endif
       if (.not.ok_kepsbc) then
         write(*,*) 'The parameter kepsbc is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_varalb) then
         write(*,*) 'The parameter varalb is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_skin) then
         write(*,*) 'The parameter skin is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_massflux) then
         write(*,*) 'The parameter massflux is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_ifrad) then
         write(*,*) 'The parameter ifrad is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_sedim) then
         write(*,*) 'The parameter sedim is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_M) then
         write(*,*) 'The parameter M is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_Mice) then
         write(*,*) 'The parameter Mice is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_ns) then
         write(*,*) 'The parameter ns is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_nstep_keps) then
         write(*,*) 'The parameter nstep_keps is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_nscreen) then
         write(*,*) 'The parameter nscreen is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_SoilType) then
         write(*,*) 'The parameter SoilType is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_Tinitlength) then
         write(*,*) 'The parameter T_profile is missing in setup file'
         all_par_present = .false.
       endif
       if (.not.ok_dyn_pgrad) then
         write(*,*) 'The parameter dyn_pgrad is missing in setup file'
         all_par_present = .false.
       endif
       if (.not.ok_monthly_out) then
         write(*,*) 'The parameter monthly_out is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_daily_out) then
         write(*,*) 'The parameter daily_out is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_hourly_out) then
         write(*,*) 'The parameter hourly_out is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_time_series) then
         write(*,*) 'The parameter time_series is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_everystep) then
         write(*,*) 'The parameter everystep is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_turb_out) then
         write(*,*) 'The parameter turb_out is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_scale_output) then
         write(*,*) 'The parameter scale_output is missing in setup file'
         all_par_present = .false.
       endif
       if (.not.ok_ngrid_out) then
         write(*,*) 'The parameter ngrid_out is missing in setup file'         
         all_par_present = .false.
       endif
       if (.not.ok_ngridsoil_out) then
         write(*,*) 'The parameter ngridsoil_out is missing in setup file'         
         all_par_present = .false.
       endif  
       if (.not.ok_assim) then
         write(*,*) 'The parameter assim is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_error_cov) then
         write(*,*) 'The parameter error_cov is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_path) then
         write(*,*) 'The parameter path is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_accum_begin) then
         write(*,*) 'The parameter accum_begin is missing in setup file'
         all_par_present = .false.
       endif         
       if (.not.ok_accum_end) then
         write(*,*) 'The parameter accum_end is missing in setup file'
         all_par_present = .false.
       endif  
       if (.not.ok_zero_model) then
         write(*,*) 'The parameter zero_model is missing in setup file'
         all_par_present = .false.
       endif      
       
       
       if (.not.all_par_present) then
         write(*,*) 'Not all necessary parameters are present in setup file: STOP'
         STOP
       endif       
       
       firstcall = .false. 
       
      endif

      END SUBROUTINE DEFINE_driving_params


      SUBROUTINE READPAR(LINE)
      implicit none

      character(len=200) :: line
      character(len=10) :: chardate
      integer(4) :: i, j
      integer(4) :: n2
      integer(4) :: n1
      integer(4) :: ncol
      real(8), allocatable :: work(:,:)
      
      SAVE

      i=1
      if (line(1:1)==' '.or.line(1:1)==char(1)) then
        do while (line(i:i)==' '.or.line(i:i)==char(1))
          i=i+1
        enddo
        n1=i
        do while (line(i:i)/=' '.and.line(i:i)/=char(1))
          i=i+1
        enddo
        n2=i-1
      else
        n1=1
        do while (line(i:i)/=' '.and.line(i:i)/=char(1))
          i=i+1
        enddo
        n2=i-1
      endif
      
      if (line(n1:n1)/='#') then

!     The group of controls for model physics

       SELECT CASE (line(n1:n2))
       CASE ('PBLpar')
        PBLpar   = igetvarval(n1,n2,line,'Param. of surf. fluxes   ')
        ok_PBLpar = .true.
       CASE ('runmode')
        runmode  = igetvarval(n1,n2,line,'Mode of model run        ')
        ok_runmode = .true.
       CASE ('relwind')
        relwind  = igetvarval(n1,n2,line,'Switch for relative wind ')
        ok_relwind = .true.
       CASE ('waveenh')
        waveenh  = igetvarval(n1,n2,line,'Switch for shal.wat.corr.')
        ok_waveenh = .true.
       CASE ('kwe')
        kwe      = getvarval (n1,n2,line,'The wave break TKE-effect')
        ok_kwe = .true.
       CASE ('Turbpar')
        Turbpar  = igetvarval(n1,n2,line,'Param. of turbulence     ')
        ok_Turbpar = .true.
       CASE ('stabfunc')
        stabfunc = igetvarval(n1,n2,line,'Stability function is    ')
        ok_stabfunc = .true.
       CASE ('kepsbc')
        kepsbc   = igetvarval(n1,n2,line,'K-EPS bound. conditions  ')
        ok_kepsbc = .true.        
       CASE ('varalb')
        varalb   = igetvarval(n1,n2,line,'Switch for variab. albedo')
        ok_varalb = .true.
       CASE ('soiltype')
        soiltype = igetvarval(n1,n2,line,'Soil type                ')
        ok_soiltype = .true.
       CASE ('soil_depth')
        depth    = getvarval (n1,n2,line,'Soil depth, m            ')
        ok_depth = .true.
       CASE ('skin')
        skin     = igetvarval(n1,n2,line,'Switch for skin at surf. ')
        ok_skin = .true.
       CASE ('sedim')
        sedim    = igetvarval(n1,n2,line,'Switch for sedimentation ')
        ok_sedim = .true.
       CASE ('massflux')
        massflux = igetvarval(n1,n2,line,'Switch for massflux conv.')
        ok_massflux = .true.
       CASE ('ifrad')
        ifrad    = igetvarval(n1,n2,line,'Switch for radiation     ')
        ok_ifrad = .true.
       CASE ('sensflux0')
        sensflux0= getvarval (n1,n2,line,'Const. sens. flux, W/m**2')
        ok_sensflux0 = .true.
       CASE ('momflux0')
        momflux0 = getvarval (n1,n2,line,'Const. mom. flux,  N/m**2')
        ok_momflux0 = .true.
       CASE ('dyn_pgrad')
        dyn_pgrad= igetvarval(n1,n2,line,'Dynamic pressure gradient')
        ok_dyn_pgrad = .true. 
       CASE ('zero_model')
        zero_model &
        &        = igetvarval(n1,n2,line,'Zero-dimensional model   ')
        ok_zero_model = .true.

       CASE ('path')
        read (line((n2+1):100),*) path
!        print*, 'path=', path
        ok_path = .true.

!     The group of numerical scheme parameters
       CASE ('nstep_keps')
        nstep_keps &
        &        = igetvarval(n1,n2,line,'Number of k-eps timesteps')
        ok_nstep_keps = .true.
       CASE ('M')
        M        = igetvarval(n1,n2,line,'Number of layers in water')
        ok_M = .true.
       CASE ('Mice')
        Mice     = igetvarval(n1,n2,line,'Number of layers in ice  ')
        ok_Mice = .true.
       CASE ('ns')
        ns       = igetvarval(n1,n2,line,'Number of layers in soil ')
        ok_ns = .true.
       CASE ('d_surf')
        d_surf   = getvarval (n1,n2,line,'Surface zooming,     n/d ')
        ok_d_surf = .true.
       CASE ('d_bot')
        d_bot    = getvarval (n1,n2,line,'Bottom zooming,      n/d ')
        ok_d_bot = .true.

!     The group of initial conditions

       CASE ('T_profile')
        Tinitlength &
        &        =igetvarval (n1,n2,line,'Init. T-profile length   ') 
        allocate (Tinitprof  (1:Tinitlength) )
        allocate (zTinitprof (1:Tinitlength) )
        ncol = 2
        allocate (work(Tinitlength,ncol))
        call READPROFILE(nunit,Tinitlength,ncol,work) 
        zTinitprof(:) = work (:,1)
        Tinitprof (:) = work (:,2)
        deallocate(work)
        ok_Tinitlength = .true.
          
!      The group of output controls

       CASE ('nscreen')
        nscreen  = igetvarval(n1,n2,line,'Screen output interval   ')
        ok_nscreen = .true.
       CASE ('monthly')
        monthly_out &
        &        = igetvarval(n1,n2,line,'Monthly output switch    ')
        ok_monthly_out = .true.
       CASE ('daily')
        daily_out= igetvarval(n1,n2,line,'Daily output switch      ')
        ok_daily_out = .true.
       CASE ('hourly')
        hourly_out &
        &        = igetvarval(n1,n2,line,'Hourly output switch     ')
        ok_hourly_out = .true.
       CASE ('everystep')
        everystep &
        &        = igetvarval(n1,n2,line,'Everystep output switch  ')
        ok_everystep = .true.
       CASE ('dt_out')
        dt_out   = getvarval (n1,n2,line,'Output interval, hours   ')
        ok_dt_out = .true.
       CASE ('time_series')
        time_series &
        &        = igetvarval(n1,n2,line,'Time series output switch')
        ok_time_series = .true.
       CASE ('turb_out')
        turb_out = igetvarval(n1,n2,line,'Switch turb. char. output')
        ok_turb_out = .true.
       CASE ('scale_output')
        scale_output &
        &        = igetvarval(n1,n2,line,'Switch turb. out. scaling')
        ok_scale_output = .true.
       CASE ('accum_begin')
        j        = igetvarval(n1,n2,line,'Accumulation period begin')
        ok_accum_begin = .true.
        write (chardate,'(i10)') j
        read (chardate,'(i4,3i2)') year_accum_begin, month_accum_begin, &
        & day_accum_begin, hour_accum_begin
       CASE ('accum_end')
        j        = igetvarval(n1,n2,line,'Accumulation period end  ')
        ok_accum_end = .true.
        write (chardate,'(i10)') j
        read (chardate,'(i4,3i2)') year_accum_end, month_accum_end, &
        & day_accum_end, hour_accum_end 
       CASE ('ngrid_out')
        ngrid_out &
        &        = igetvarval(n1,n2,line,'Switch for output regrid ')
        allocate (zgrid_out  (1:ngrid_out) )
        ok_ngrid_out = .true.
        if (ngrid_out > 0) then
          ncol = 1
          allocate (work(1:ngrid_out,ncol))
          call READPROFILE(nunit,ngrid_out,ncol,work) 
          zgrid_out(:) = work (:,1)
          deallocate(work)
        endif
       CASE ('ngridsoil_out')
        ngridsoil_out &
        &        = igetvarval(n1,n2,line,'Switch soil out. regrid  ')
        allocate (zgridsoil_out  (1:ngridsoil_out) )
        ok_ngridsoil_out = .true.
        if (ngridsoil_out > 0) then
          ncol = 1
          allocate (work(1:ngridsoil_out,ncol))
          call READPROFILE(nunit,ngridsoil_out,ncol,work) 
          zgridsoil_out(:) = work (:,1)
          deallocate(work)
        endif

!     The group of physical properties

!     The group of data assimilation controls

       CASE ('assim')
        assim    = igetvarval(n1,n2,line,'The assimilation switch  ')
        ok_assim = .true.
       CASE ('as_window')
        as_window= igetvarval(n1,n2,line,'Assimilation window      ')
        ok_as_window = .true.
       CASE ('error_cov')
        error_cov= igetvarval(n1,n2,line,'The switch for err_cov   ')
        ok_error_cov = .true.

!     The group of parameters of inflows and outflows

       CASE ('tribheat')
        tribheat = igetvarval(n1,n2,line,'Switch for trib. heat    ')
        ok_tribheat = .true.
       CASE ('inflowprof')
          N_tribin = 1
          allocate (width_tribin  (1:N_tribin ,1:M+1) )
          allocate (U_tribin      (1:N_tribin ,1:M+1) )
          allocate (T_tribin      (1:N_tribin ,1:M+1) )
          ncol = 4
          allocate (work(M+1,ncol))
          call READPROFILE(nunit,M+1,ncol,work) 
!         In current version there is only ONE inflow     
          width_tribin(1,:) = work (:,2)
          U_tribin    (1,:) = work (:,3)
          T_tribin    (1,:) = work (:,4)
          deallocate(work)
!          ok_inflowprof = .true.
       CASE ('outflowprof')
          N_tribout = 1
          allocate (width_tribout  (1:N_tribout ,1:M+1) )
          allocate (U_tribout      (1:N_tribout ,1:M+1) )
          ncol = 3
          allocate (work(M+1,ncol))
          print*, 'The profiles in inflow are done'
          call READPROFILE(nunit,M+1,ncol,work) 
!         In current version there is only ONE outflow 
          width_tribout(1,:) = work (:,2)
          U_tribout    (1,:) = work (:,3)
          print*, 'The profiles in ouflow are done'
          deallocate(work)
!          ok_outflowprof = .true.
        
       CASE DEFAULT
        if (.not.line(n1:n2)=='end') then
         print*,'Unknown keyword [',line(n1:n2),'] in setup file: STOP'
         STOP
        endif
       END SELECT

      endif

      END SUBROUTINE READPAR


      SUBROUTINE READPROFILE(iunit,N_rows,N_coloumns,work)
      implicit none

!     Input variables:
      integer(4), intent(in) :: iunit
      integer(4), intent(in) :: N_rows
      integer(4), intent(in) :: N_coloumns

!     Output variables:
      real(8), intent(out) :: work(N_rows, N_coloumns)

!     Local variables
      integer(4) :: i
      integer(4) :: j

      do i=1, N_rows
        read(iunit,*) (work(i,j), j=1,N_coloumns)
      enddo

      END SUBROUTINE READPROFILE


      END MODULE DRIVING_PARAMS




      MODULE ARRAYS
      
!     MODULE ARRAYS contains allocatable arrays of the model  

      use DRIVING_PARAMS
      
      integer(4), parameter :: soil_indic = 1
      integer(4), parameter :: ice_indic = 2
      integer(4), parameter :: snow_indic = 3
      integer(4), parameter :: water_indic = 4
      integer(4), parameter :: deepice_indic = 5
      integer(4), parameter :: water_salinity_indic = 6
      integer(4), parameter :: soil_salinity_indic = 7
      integer(4), parameter :: water_methane_indic = 8
      integer(4), parameter :: water_oxygen_indic = 9
      integer(4), parameter :: water_carbdi_indic = 10

      integer(4) snow,ice,water,deepice,nstep 
      real(8) time,tempairf,relhumf,humairf,pressuref, &
      & windf,shortwavef,longwavef,precipf,time_old,tempair_old, &
      & relhum_old,humair_old,pressure_old,wind_old, &
      & shortwave_old,longwave_old,precip_old,Erad, &
      & wdir_old,wdirf,dep_av
      real(8) zinv
      integer(4) time_int,time_int_old,man,par,flag_snow_init

!     Bathymetry group
      real(8), allocatable :: area_int(:), area_half(:)
 
!     Methane group
      real(8) :: veg
      real(8) :: fplant, fdiff, fdiff_lake_surf
      real(8) :: febul
!      real(8) :: febul2(1:2) ! two-meth
!      real(8) :: fdiff2(1:2), fdiff_lake_surf2(1:2) ! two-meth
      real(8) :: febultot(1:2), fdifftot(1:2), fdiff_lake_surftot(1:2)
!      real(8) :: febultot2(1:2,1:2), fdifftot2(1:2,1:2), fdiff_lake_surftot2(1:2,1:2) ! two-meth
      real(8) :: plant_sum, bull_sum, oxid_sum, rprod_sum
      real(8) :: rprod_total_oldC, rprod_total_newC
      real(8) :: rprod_total_newC_integr(1:2), rprod_total_oldC_integr(1:2)
!      real(8) :: rprod_total_newC_integr2(1:2,1:2), rprod_total_oldC_integr2(1:2,1:2) ! two-meth
      real(8) :: h_talik
      real(8) :: anox
      real(8), allocatable :: rootss(:), qmethane(:), TgrAnn(:)
      real(8), allocatable :: qsoil(:), qwater(:,:)
!      real(8), allocatable :: qsoil2(:,:), qwater2(:,:,:)  ! two-meth
      real(8), allocatable :: lammeth(:)
      real(8) :: tot_ice_meth_bubbles
!      real(8) :: tot_ice_meth_bubbles2(1:2)  ! two-meth
      
!     Oxygen group      
      real(8), allocatable :: oxyg(:,:)
      real(8), allocatable :: lamoxyg(:)
      
!     Carbon dioxide group 
      real(8), allocatable :: carbdi(:,:)
      real(8), allocatable :: lamcarbdi(:)
      
      real(8) totalevaps,totalmelts,totalprecips
      real(8), allocatable :: Tskin(:)
      real(8) :: T_0dim
      real(8) h1,l1,hs1,ls1
      real(8) dhw,dhw0,dhi,dhi0,dls,dls0,dhwfsoil
      real(8) :: Buoyancy0,H_mixed_layer,w_conv_scale,T_conv_scale,H_entrainment
      real(8) :: S_integr_positive, S_integr_negative, Gen_integr, &
      & eps_integr, E_integr
      real(8) :: Seps_integr_positive, Seps_integr_negative, &
      & Geneps_integr, epseps_integr
      real(8) :: TKE_balance, eps_balance
      real(8) :: SR_botsnow
      real(8) :: dt_keps
      real(8) :: dt05, dt_inv, dt_inv05
      real(8) :: deltaskin
      real(8), allocatable:: dz_full(:), z_full(:), z_half(:)
      real(8), allocatable :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
      real(8), allocatable :: ddzi(:), ddzi05(:)
      real(8), allocatable :: dzeta_int(:), dzeta_05int(:)
      real(8), allocatable :: dzetai_int(:), dzetai_05int(:)      
      real(8), allocatable :: Tw1(:),Tw2(:),Ti1(:),Ti2(:),Tis1(:), &
      & Tis2(:),RS(:),SR(:),lamw(:),lamsal(:),Sal1(:),Sal2(:), &
      & Sals1(:),Sals2(:),SRi(:),SRdi(:)
      real(8), allocatable:: Tsoil1(:),rosoil(:),csoil(:), &
      & lamsoil(:), dzs(:), dzss(:), wi1(:), wl1(:), Tsoil2(:), &
      & zsoil(:), wsoil(:) 
      real(8), allocatable :: z_watersoil(:)
      real(8), allocatable:: S(:),Gen(:),F(:),TKE_turb_trans(:)
      real(8), allocatable:: KT(:),k2(:),u1(:),v1(:),w1(:),E1(:),eps1(:)
      real(8), allocatable:: WLM0(:),WLM7(:),bH(:),PSIMAX(:),POR(:),FLWMAX(:),DLMAX(:)
      real(8), allocatable:: KC(:),KLengT(:),RSR(:) 
      real(8), allocatable:: dep_2d(:,:)
      integer(4), allocatable :: init(:,:), num(:)
      real(8), allocatable:: Tsoil4(:), Tsoil3(:), &
      & wl2(:), wl3(:), wl4(:), wl5(:),  &
      & wi2(:), wi3(:), lammoist(:), rosdry(:), filtr(:), &
      & wa(:)
      real(8), allocatable::E2(:),eps2(:),eps3(:), &
      & u2(:),v2(:),w2(:),k1(:),k3(:),k4(:),k5(:),u3(:), &
      & v3(:),C1aup(:),Re(:),row(:),Ri(:),uv(:),E_it1(:), &
      & E_it2(:),C_num(:),E_it3(:),Feps(:),E_it21(:), &
      & eps_it21(:),eps_it1(:),eps_it2(:), &
      & l(:),k2_mid(:),k3_mid(:),k4_mid(:),k5_mid(:), &
      & Um(:),Vm(:),E12(:),eps12(:),knum(:),k2t(:), &
      & Eeps(:),Ecent(:),Epscent(:),res_E(:), &
      & res_eps(:),dresd(:,:,:,:),dres2dE(:),dres2deps(:), &
      & veg_sw(:) 
      real(8), allocatable:: WU_(:),WV_(:),GAMT(:),GAMU(:),GAMV(:), TF(:),KLengu(:)
      real(8), allocatable:: PEMF    (:) , PDENS_DOWN (:), &
      &                       PT_DOWN (:) , PSAL_DOWN  (:), &
      &                       pt_down_f(:)
      real(8), allocatable:: k_turb_T_flux(:), T_massflux(:)

      data nstep /0/

      SAVE

      END MODULE ARRAYS

      

      MODULE TURB_CONST

!     PHYSICAL CONSTANTS
       real(8), parameter:: kar       = 0.4d0
       real(8), parameter:: niu       = 1.007d-6

!     DIMENSIONLESS CONSTANTS IN E-EPS PARAMETERIZATION, according to Goudsmit et al., 2002
       real(8), parameter:: CE0       = 0.09d0 
       real(8), parameter:: CEt0      = 0.072d0
       real(8), parameter:: ceps1     = 1.44d0
       real(8), parameter:: ceps2     = 1.92d0
       real(8), parameter:: sigmaE    = 1.d0
       real(8), parameter:: sigmaeps  = 1.3d0
!       real(8), parameter:: sigmaeps  = 1.111d0 ! according to Burchard, 2002

       real(8), parameter:: CL = CE0**(0.75d0) ! This value results from assumption
                                               ! that shear production is balanced by dissipation
                                               ! (Burchard, 2002)
                                               
       real(8), parameter :: C_wstar = 0.4d0   ! Deardorff, 1970

       real(8), parameter:: lam_T     = 1.d0 
       real(8), parameter:: lam_gen   = 1.d0 
       real(8), parameter:: lamTM     = 0.14d0

!	   real(8), parameter:: cmu_0     = 0.094d0
!      Other values for cmu_0, proposed in literature, are
!      0.121, 0.077 - see Burchard, 2002

       real(8), parameter:: lam_E0     = CE0/sigmaE   
       real(8), parameter:: lam_eps0   = CE0/sigmaeps 

!      The coefficients in formulas for stability functions (Canuto et al., 2001)
       real(8), parameter:: CEcoef01   = 0.127d0
       real(8), parameter:: CEcoef02   = 0.1526d-1
       real(8), parameter:: CEcoef03   = 0.16d-3

       real(8), parameter:: CEtcoef01  = 0.1190d0
       real(8), parameter:: CEtcoef02  = 0.4294d-2
       real(8), parameter:: CEtcoef03  = 0.66d-3

       real(8), parameter:: CEcoef11   = 1.d0
       real(8), parameter:: CEcoef12   = 0.2d0
       real(8), parameter:: CEcoef13   = 0.315d-1
       real(8), parameter:: CEcoef14   = 0.58d-2
       real(8), parameter:: CEcoef15   = 0.4d-2
       real(8), parameter:: CEcoef16   = 0.4d-4
 
!      The parameters of Mellor and Yamada model (1982)	 

       real(8), parameter:: A1  =  0.92d0
       real(8), parameter:: B1  =  1.66d+1
       real(8), parameter:: C1  =  0.8d-1
       real(8), parameter:: A2  =  0.74d0
       real(8), parameter:: B2  =  1.01d+1
       real(8), parameter:: CL_K_KL_MODEL = 2.d0**(1.5d0)/B1

!     Co is according to Satyanarayana et al., 1999
       real(8), parameter:: Co        = 1.9d0

!     CONSTANTS IN CONVECTION PARAMETERIZATION
       real(8), parameter:: KC0       = 0.d0 !1.d+5 !500.
       real(8), parameter:: dtdz0     = 2.d0

!     CONSTANTS OF SIMOES'S PARAMETERIZATION
       real(8), parameter:: Cs        = 1.5d0  !0.15
       real(8), parameter:: Cs1       = 100.d0    
      
!     lo=0.4*ddz*h1*3.                                                         VS,06.2007
       real(8), parameter:: L0        = 0.1d0 !0.1d0

       real(8), parameter:: CON0      = 0.05d0
       real(8), parameter:: CON1      = 0.11d0
       real(8), parameter:: CON2      = 1.56d0
       real(8), parameter:: CONUV     = 1.d0

      END MODULE TURB_CONST
      
      
!      MODULE ASSIM_VAR
!      use DRIVING_PARAMS
!      integer n_modvar, n_obsvar, size_mod, size_obs, ntim
!      
!      real(8), allocatable:: y_obs(:,:), x_mod(:,:), b_mod(:,:), &
!      & r_obs(:,:), h_obs(:,:), covar(:,:,:,:), c_diff(:,:,:), &
!      & RpHBHT(:,:), dep_obs(:,:),y_obs_mean(:),x_mod_mean(:), &
!      & delta(:)
!
!      SAVE
!     
!      contains
!      subroutine alloc_assim_var
!      n_modvar = 1
!      n_obsvar = 1
!      size_mod = (M+1)*n_modvar
!      size_obs = (M+1)*n_obsvar
!      
!      allocate (y_obs(1:ntim, 1:size_obs), &
!      &          x_mod(1:ntim, 1:size_mod), &
!      &          dep_obs(1:ntim, 1:size_obs))
!      allocate (r_obs(1:size_obs, 1:size_obs), &
!      &          RpHBHT(1:size_obs,1:size_obs), &
!      &          b_mod(1:size_mod, 1:size_mod), &
!      &          h_obs(1:size_obs, 1:size_mod))
!      allocate (c_diff(1:M+1,1:n_obsvar,1:n_obsvar))
!      allocate (covar(1:M+1, 1:M+1, 1:n_obsvar, 1:n_obsvar))
!      allocate (y_obs_mean(1:size_obs), &
!      &          x_mod_mean(1:size_mod), &
!      &          delta     (1:size_obs))
!      
!      y_obs=0.
!      x_mod=0.
!      r_obs=0.
!      b_mod=0. 
!      h_obs=0.
!      c_diff=0.
!      covar=0.
!      end subroutine alloc_assim_var
!      
!      END MODULE ASSIM_VAR


MODULE EVOLUTION_VARIABLES

real(8), allocatable, private :: l1_2d(:,:)
real(8), allocatable, private :: h1_2d(:,:)
real(8), allocatable, private :: ls1_2d(:,:)
real(8), allocatable, private :: hs1_2d(:,:)
real(8), allocatable, private :: time_2d(:,:) 
real(8), allocatable, private :: cdm_2d(:,:)
real(8), allocatable, private :: u_2d(:,:,:)
real(8), allocatable, private :: v_2d(:,:,:)
real(8), allocatable, private :: Tsoil1_2d(:,:,:)
real(8), allocatable, private :: wi1_2d(:,:,:)
real(8), allocatable, private :: wl1_2d(:,:,:)
real(8), allocatable, private :: Tw1_2d(:,:,:)
real(8), allocatable, private :: Ti1_2d(:,:,:)
real(8), allocatable, private :: Tis1_2d(:,:,:)
real(8), allocatable, private :: dz_2d(:,:,:)
real(8), allocatable, private :: T_2d(:,:,:)
real(8), allocatable, private :: wl_2d(:,:,:)
real(8), allocatable, private :: dens_2d(:,:,:)
real(8), allocatable, private :: E_2d(:,:,:)
real(8), allocatable, private :: eps_2d(:,:,:)
real(8), allocatable, private :: dhwfsoil_2d(:,:) 
real(8), allocatable, private :: Sal1_2d(:,:,:)
real(8), allocatable, private :: Sals1_2d(:,:,:)
real(8), allocatable, private :: Elatent_2d(:,:)
real(8), allocatable, private :: dhw_2d(:,:)
real(8), allocatable, private :: dhw0_2d(:,:)
real(8), allocatable, private :: dhi_2d(:,:)
real(8), allocatable, private :: dhi0_2d(:,:)
real(8), allocatable, private :: dls0_2d(:,:)
real(8), allocatable, private :: velfrict_2d(:,:)
real(8), allocatable, private :: roughness_2d(:,:)
real(8), allocatable, private :: eflux0_kinem_2d(:,:)
real(8), allocatable, private :: lamw_2d(:,:,:)
real(8), allocatable, private :: snmelt_2d(:,:)
real(8), allocatable, private :: Tskin_2d(:,:)
real(8), allocatable, private :: k_turb_T_flux_2d(:,:,:)
real(8), allocatable, private :: qwater_2d(:,:,:)
real(8), allocatable, private :: qsoil_2d(:,:,:)
real(8), allocatable, private :: oxyg_2d(:,:,:)
real(8), allocatable, private :: carbdi_2d(:,:,:)
real(8), allocatable, private :: tot_ice_meth_bubbles_2d(:,:)

integer(4), allocatable, private :: fl_sn_2d(:,:)
integer(4), allocatable, private :: fl_sn_init_2d(:,:) 
integer(4), allocatable, private :: itop_2d(:,:)
integer(4), allocatable, private :: nstep_2d(:,:)

logical, private :: arrays_allocated = .false.

SAVE

contains
SUBROUTINE ALLOCATE_ARRAYS(nx, ny, M, Mice, ns, ms, ml)

implicit none

! Input variables
 integer(4), intent(in) :: nx, ny
 integer(4), intent(in) :: M, Mice, ns, ms, ml

 allocate ( l1_2d(nx,ny) ) ; l1_2d = 0
 allocate ( h1_2d(nx,ny) ) ; h1_2d = 0
 allocate ( ls1_2d(nx,ny) ) ; ls1_2d = 0
 allocate ( cdm_2d(nx,ny) ) ; cdm_2d = 0
 allocate ( u_2d(M+1,nx,ny) ) ; u_2d = 0
 allocate ( v_2d(M+1,nx,ny) ) ; v_2d = 0
 allocate ( E_2d(M+1,nx,ny) ) ; E_2d = 0
 allocate ( eps_2d(M+1,nx,ny) ) ; eps_2d = 0
 allocate ( Tsoil1_2d(ns+1,nx,ny) ) ; Tsoil1_2d = 0
 allocate ( wi1_2d(ns+1,nx,ny) ) ; wi1_2d = 0
 allocate ( wl1_2d(ns+1,nx,ny) ) ; wl1_2d = 0
 allocate ( Tw1_2d(M+1,nx,ny) ) ; Tw1_2d = 0
 allocate ( Sal1_2d(M+1,nx,ny) ) ; Sal1_2d = 0
 allocate ( Ti1_2d(Mice+1,nx,ny) ) ; Ti1_2d = 0
 allocate ( Tis1_2d(Mice+1,nx,ny) ) ; Tis1_2d = 0
 allocate ( k_turb_T_flux_2d(M+1,nx,ny) ) ; k_turb_T_flux_2d = 0
 allocate ( Sals1_2d(ns,nx,ny) ) ; Sals1_2d = 0
 allocate ( hs1_2d(nx,ny) ) ; hs1_2d = 0
 allocate ( dz_2d(ms,nx,ny) ) ; dz_2d = 0
 allocate ( T_2d(ml,nx,ny) ) ; T_2d = 0
 allocate ( wl_2d(ml,nx,ny) ) ; wl_2d = 0
 allocate ( dens_2d(ms,nx,ny) ) ; dens_2d = 0
 allocate ( time_2d(nx,ny) ) ; time_2d = 0
 allocate ( dhwfsoil_2d(nx,ny) ) ; dhwfsoil_2d = 0
 allocate ( Elatent_2d(nx,ny) ) ; Elatent_2d = 0
 allocate ( dhw_2d(nx,ny) ) ; dhw_2d = 0
 allocate ( dhw0_2d(nx,ny) ) ; dhw0_2d = 0
 allocate ( dhi_2d(nx,ny) ) ; dhi_2d = 0
 allocate ( dhi0_2d(nx,ny) ) ; dhi0_2d = 0
 allocate ( dls0_2d(nx,ny) ) ; dls0_2d = 0
 allocate ( velfrict_2d(nx,ny) ) ; velfrict_2d = 0
 allocate ( roughness_2d(nx,ny) ) ; roughness_2d = 0
 allocate ( eflux0_kinem_2d(nx,ny) ) ; eflux0_kinem_2d = 0
 allocate ( lamw_2d(nx,ny,M) ) ; lamw_2d = 0
 allocate ( snmelt_2d(nx,ny) ) ; snmelt_2d = 0
 allocate ( fl_sn_2d(nx,ny) ) ; fl_sn_2d = 0
 allocate ( fl_sn_init_2d(nx,ny) ) ; fl_sn_init_2d = 0
 allocate ( itop_2d(nx,ny) ) ; itop_2d = 0
 allocate ( nstep_2d(nx,ny) ) ; nstep_2d = 0
 allocate ( Tskin_2d(nx,ny) ) ; Tskin_2d = 0
 allocate ( qwater_2d(M+1,nx,ny) ) ; qwater_2d = 0.
 allocate ( qsoil_2d(ns,nx,ny) ) ; qsoil_2d = 0.
 allocate ( oxyg_2d(M+1,nx,ny) ) ; oxyg_2d = 0.
 allocate ( carbdi_2d(M+1,nx,ny) ) ; carbdi_2d = 0. 
 allocate ( tot_ice_meth_bubbles_2d(nx,ny) ) ; tot_ice_meth_bubbles_2d = 0.

 arrays_allocated = .true.

END SUBROUTINE ALLOCATE_ARRAYS



SUBROUTINE UPDATE_CURRENT_TIMESTEP ( &
& ix, iy, nx, ny, M, Mice, ns, ms, ml, &
& l1, h1, ls1, hs1, &
& u1, v1, &
& E1, eps1, k_turb_T_flux, &
& Tsoil1, Sals1, wi1, wl1, &
& Tw1, Sal1, lamw, &
& Tskin, &
& Ti1, Tis1, &
& dz, T, wl, dens, &
& qwater, qsoil, &
& oxyg, carbdi, &
& snmelt, &
& cdmw, &
& time, &
& dhwfsoil, &
& Elatent, &
& dhw, &
& dhw0, &
& dhi, &
& dhi0, &
& dls0, &
& velfrict, &
& roughness, &
& eflux0_kinem, &
& tot_ice_meth_bubbles, &

& flag_snow, flag_snow_init, &
& itop, &
& nstep )

implicit none

! Input variables
integer(4), intent(in) :: ix, iy, nx, ny
integer(4), intent(in) :: M, Mice, ns, ms, ml

! Output variables
real(8), intent(out) :: l1, h1, ls1, hs1
real(8), intent(out) :: u1(1:M+1), v1(1:M+1)
real(8), intent(out) :: E1(1:M+1), eps1(1:M+1), k_turb_T_flux(1:M+1)
real(8), intent(out) :: Tsoil1(1:ns), Sals1(1:ns), wi1(1:ns), wl1(1:ns)
real(8), intent(out) :: Tw1(1:M+1), Sal1(1:M+1), lamw(1:M)
real(8), intent(out) :: Tskin
real(8), intent(out) :: Ti1(1:Mice+1), Tis1(1:Mice+1)
real(8), intent(out) :: dz(1:ms), T(1:ml), wl(1:ml), dens(1:ms)
real(8), intent(out) :: qwater(1:M+1), qsoil(1:ns)
real(8), intent(out) :: oxyg(1:M+1), carbdi(1:M+1)
real(8), intent(out) :: snmelt
real(8), intent(out) :: cdmw
real(8), intent(out) :: time
real(8), intent(out) :: dhwfsoil
real(8), intent(out) :: Elatent
real(8), intent(out) :: dhw
real(8), intent(out) :: dhw0
real(8), intent(out) :: dhi
real(8), intent(out) :: dhi0
real(8), intent(out) :: dls0
real(8), intent(out) :: velfrict
real(8), intent(out) :: roughness
real(8), intent(out) :: eflux0_kinem
real(8), intent(out) :: tot_ice_meth_bubbles

integer(4), intent(out) :: flag_snow, flag_snow_init
integer(4), intent(out) :: itop
integer(4), intent(out) :: nstep

! Local variables
integer(4) :: i

if (.not.arrays_allocated) &
& call ALLOCATE_ARRAYS(nx, ny, M, Mice, ns, ms, ml)

l1 = l1_2d(ix,iy) 
h1 = h1_2d(ix,iy)
ls1 = ls1_2d(ix,iy) 
hs1 = hs1_2d(ix,iy)

do i = 1, M+1
  u1(i) = u_2d(i,ix,iy)
  v1(i) = v_2d(i,ix,iy)
enddo

do i = 1, M
  E1(i) = E_2d(i,ix,iy)
  eps1(i) = eps_2d(i,ix,iy)
  k_turb_T_flux(i) = k_turb_T_flux_2d(i,ix,iy)
enddo
    
do i = 1, ns 
  Tsoil1(i) = Tsoil1_2d(i,ix,iy)
  Sals1(i) = Sals1_2d(i,ix,iy)
  wi1(i) = wi1_2d(i,ix,iy)
  wl1(i) = wl1_2d(i,ix,iy)
  qsoil(i) = qsoil_2d(i,ix,iy)
enddo
    
do i = 1, M+1 
  Tw1(i) = Tw1_2d(i,ix,iy)
  Sal1(i) = Sal1_2d(i,ix,iy)
  qwater(i) = qwater_2d(i,ix,iy)
  oxyg(i) = oxyg_2d(i,ix,iy)
  carbdi(i) = carbdi_2d(i,ix,iy)
enddo
  
Tskin = Tskin_2d(ix,iy)
 
do i = 1, Mice + 1 
  Ti1(i) = Ti1_2d(i,ix,iy)
  Tis1(i) = Tis1_2d(i,ix,iy)
enddo

flag_snow = fl_sn_2d(ix,iy)
flag_snow_init = fl_sn_init_2d(ix,iy)

itop = itop_2d(ix,iy)

do i = max(1,itop), ms
  dz(i) = dz_2d(i,ix,iy)
  T(i) = T_2d(i,ix,iy)
  wl(i) = wl_2d(i,ix,iy)
  dens(i) = dens_2d(i,ix,iy)
enddo

snmelt = snmelt_2d(ix,iy)
 
cdmw = cdm_2d(ix,iy)
  
time = time_2d(ix,iy)
nstep = nstep_2d(ix,iy)
dhwfsoil = dhwfsoil_2d(ix,iy)
Elatent = Elatent_2d(ix,iy)
dhw = dhw_2d(ix,iy)
dhw0 = dhw0_2d(ix,iy)
dhi = dhi_2d(ix,iy)
dhi0 = dhi0_2d(ix,iy)
dls0 = dls0_2d(ix,iy)
velfrict = velfrict_2d(ix,iy)
roughness = roughness_2d(ix,iy)
eflux0_kinem = eflux0_kinem_2d(ix,iy)

tot_ice_meth_bubbles = tot_ice_meth_bubbles_2d(ix,iy)

lamw(1:M) = lamw_2d(ix,iy,1:M)
  
END SUBROUTINE UPDATE_CURRENT_TIMESTEP



SUBROUTINE UPDATE_NEXT_TIMESTEP ( &
& ix, iy, nx, ny, M, Mice, ns, ms, ml, &
& l1, h1, ls1, hs1, &
& u1, v1, &
& E1, eps1, k_turb_T_flux, &
& Tsoil1, Sals1, wi1, wl1, &
& Tw1, Sal1, lamw, &
& Tskin, &
& Ti1, Tis1, &
& dz, T, wl, dens, &
& qwater, qsoil, &
& oxyg, carbdi, &
& snmelt, &
& cdmw, &
& time, &
& dhwfsoil, &
& Elatent, &
& dhw, &
& dhw0, &
& dhi, &
& dhi0, &
& dls0, &
& velfrict, &
& roughness, &
& eflux0_kinem, &
& tot_ice_meth_bubbles, &

& flag_snow, flag_snow_init, &
& itop, &
& nstep )

implicit none

! Input variables
integer(4), intent(in) :: ix, iy, nx, ny
integer(4), intent(in) :: M, Mice, ns, ms, ml

real(8), intent(in) :: l1, h1, ls1, hs1
real(8), intent(in) :: u1(1:M+1), v1(1:M+1)
real(8), intent(in) :: E1(1:M+1), eps1(1:M+1), k_turb_T_flux(1:M+1)
real(8), intent(in) :: Tsoil1(1:ns), Sals1(1:ns), wi1(1:ns), wl1(1:ns)
real(8), intent(in) :: Tw1(1:M+1), Sal1(1:M+1), lamw(1:M)
real(8), intent(in) :: Tskin
real(8), intent(in) :: Ti1(1:Mice+1), Tis1(1:Mice+1)
real(8), intent(in) :: dz(1:ms), T(1:ml), wl(1:ml), dens(1:ms)
real(8), intent(in) :: qwater(1:M+1), qsoil(1:ns)
real(8), intent(in) :: oxyg(1:M+1), carbdi(1:M+1)
real(8), intent(in) :: snmelt
real(8), intent(in) :: cdmw
real(8), intent(in) :: time
real(8), intent(in) :: dhwfsoil
real(8), intent(in) :: Elatent
real(8), intent(in) :: dhw
real(8), intent(in) :: dhw0
real(8), intent(in) :: dhi
real(8), intent(in) :: dhi0
real(8), intent(in) :: dls0
real(8), intent(in) :: velfrict
real(8), intent(in) :: roughness
real(8), intent(in) :: eflux0_kinem
real(8), intent(in) :: tot_ice_meth_bubbles

integer(4), intent(in) :: flag_snow, flag_snow_init
integer(4), intent(in) :: itop
integer(4), intent(in) :: nstep

! Local variables
integer(4) :: i

if (.not.arrays_allocated) &
& call ALLOCATE_ARRAYS(nx, ny, M, Mice, ns, ms, ml)

l1_2d(ix,iy)  = l1 
h1_2d(ix,iy)  = h1
ls1_2d(ix,iy) = ls1
hs1_2d(ix,iy) = hs1

do i = 1, M+1
  u_2d(i,ix,iy) = u1(i)
  v_2d(i,ix,iy) = v1(i)
enddo

do i = 1, M
  E_2d(i,ix,iy) = E1(i)
  eps_2d(i,ix,iy) = eps1(i)
  k_turb_T_flux_2d(i,ix,iy) = k_turb_T_flux(i)
enddo
    
do i = 1, ns 
  Tsoil1_2d(i,ix,iy) = Tsoil1(i)
  Sals1_2d(i,ix,iy) = Sals1(i) 
  wi1_2d(i,ix,iy) = wi1(i)
  wl1_2d(i,ix,iy) = wl1(i)
  qsoil_2d(i,ix,iy) = qsoil(i)
enddo
    
do i = 1, M+1 
  Tw1_2d(i,ix,iy) = Tw1(i)
  Sal1_2d(i,ix,iy) = Sal1(i) 
  qwater_2d(i,ix,iy) = qwater(i)
  oxyg_2d(i,ix,iy) = oxyg(i)
  carbdi_2d(i,ix,iy) = carbdi(i)
enddo

Tskin_2d(ix,iy) = Tskin

do i = 1, Mice+1
  Ti1_2d(i,ix,iy) = Ti1(i)
  Tis1_2d(i,ix,iy) = Tis1(i)
enddo
 
fl_sn_2d(ix,iy) = flag_snow
fl_sn_init_2d(ix,iy) = flag_snow_init

itop_2d(ix,iy) = itop

do i = max(1,itop), ms
  dz_2d(i,ix,iy) = dz(i)
  T_2d(i,ix,iy) = T(i)
  wl_2d(i,ix,iy) = wl(i)
  dens_2d(i,ix,iy) = dens(i)
enddo

snmelt_2d(ix,iy) = snmelt

cdm_2d(ix,iy) = cdmw

time_2d(ix,iy) = time
nstep_2d(ix,iy) = nstep
dhwfsoil_2d(ix,iy) = dhwfsoil
Elatent_2d(ix,iy) = Elatent
dhw_2d(ix,iy) = dhw
dhw0_2d(ix,iy) = dhw0
dhi_2d(ix,iy) = dhi
dhi0_2d(ix,iy) = dhi0
dls0_2d(ix,iy) = dls0
velfrict_2d(ix,iy) = velfrict
roughness_2d(ix,iy) = roughness
eflux0_kinem_2d(ix,iy) = eflux0_kinem

tot_ice_meth_bubbles_2d(ix,iy) = tot_ice_meth_bubbles

lamw_2d(ix,iy,1:M) = lamw(1:M)

END SUBROUTINE UPDATE_NEXT_TIMESTEP

END MODULE EVOLUTION_VARIABLES


MODULE METH_OXYG_CONSTANTS

real(8), parameter :: tortuosity_coef = 0.66 ! the tortuosity of soil pores

real(8), parameter :: diff_unsat = 2.d-7 ! Diffusion constant in unsaturated soil, m**2/s
real(8), parameter :: diff_water = 2.d-9 ! Diffusion constant in water, m**2/s

real(8), parameter :: ch4_atm0 = 2.28d-6 ! concentration of methane in water, in equilibrium with
                                         ! atmospheric partial methane pressure (0.175 Pa), 
                                         ! according to Henry law with reference constant, mol/m**3
                                         ! 7.6d-5 - atmospheric concentration of methane, mol/m**3
real(8), parameter :: ch4_pres_atm0 = 1.75d-1 ! atmospheric partial methane pressure, Pa
real(8), parameter :: o2_atm0 = 2.73d-1  ! concentration of oxygen in water, in equilibrium with
                                         ! atmospheric partial oxygen pressure (21278 Pa), 
                                         ! according to Henry law with reference constant, mol/m**3
                                         ! 8. - atmospheric concentration of oxygen, mol/m**3
real(8), parameter :: o2_pres_atm0 = 21278. ! atmospheric partial oxygen pressure, Pa
real(8), parameter :: n2_atm0 = 4.96d-1  ! concentration of nitrogen in water, in equilibrium with
                                         ! atmospheric partial nitrogen pressure (80046 Pa), 
                                         ! according to Henry law with reference constant, mol/m**3
real(8), parameter :: n2_pres_atm0 = 80046. ! atmospheric partial nitrogen pressure, Pa
real(8), parameter :: co2_pres_atm0 = 39.01 ! atmospheric partial carbon dioxide pressure, Pa
real(8), parameter :: co2_atm0 = 1.3d-2  ! concentration of carbon dioxide in water, in equilibrium with
                                         ! atmospheric partial carbon dioxide pressure (39.01 Pa), 
                                         ! according to Henry law with reference constant, mol/m**3
                                        
real(8), parameter :: Henry_const0_ch4 = 1.3d-5 ! Henry constant for methane at the reference temperature, mol/(m**3*Pa)
real(8), parameter :: Henry_const0_o2 = 1.3d-5 ! Henry constant for oxygen at the reference temperature, mol/(m**3*Pa)
real(8), parameter :: Henry_const0_n2 = 6.2d-6 ! Henry constant for nitrogen at the reference temperature, mol/(m**3*Pa)
real(8), parameter :: Henry_const0_co2 = 3.36d-4 ! Henry constant for carbon doixide at the reference temperature, mol/(m**3*Pa)
real(8), parameter :: Henry_temp_ref = 298.15 ! The reference temperature for Henry constants, K
real(8), parameter :: Henry_temp_dep_ch4 = 1.7d+3 ! The temperature dependence of Henry constant for methane, K
real(8), parameter :: Henry_temp_dep_o2 = 1.7d+3 ! The temperature dependence of Henry constant for oxygen, K
real(8), parameter :: Henry_temp_dep_n2 = 1.3d+3 ! The temperature dependence of Henry constant for nitrogen, K
real(8), parameter :: Henry_temp_dep_co2 = 2.4d+3 ! The temperature dependence of Henry constant for carbon dioxide, K

real(8), parameter :: molmass_ch4 = 16. ! methane molar mass, g/mol        

real(8), parameter :: n2_exp_decay = 50. ! the decay rate in exponential law for nitrogen conecntration in soil, m**(-1)
                                         ! after (Bazhin, 2001)


!*               parameters for production
!*    =============================================================
!      real rnpp,rnppmax,rnroot,q100,r0,
!     *     vmax,rkm,q10,
!     *     rke,cmin,punveg,cthresh,
!     *     rkp,tveg,pox,rlmin,rl,rlmax
!*    =============================================================

real(8), parameter :: rnpp = 100.          ! 100. - previous value ! NPP (gC/(m**2*month))
real(8), parameter :: rnppmax = 180.       ! the maximum value of the NPP
real(8), parameter :: q100 = 6.            !
real(8), parameter :: r0 = 1.67d-7         ! the constant rate factor, mol/(m**3*s)
real(8), parameter :: r0_oliglake = 1.7d-8  !1.d-8 ! the constant rate factor for oligotrophic lake, mol/(m**3*s)
real(8), parameter :: forg0 = 0.857        ! the constant in depth dependence of methane production rate
real(8), parameter :: lambda_new_org = 3. !5. in Walter & Heimann model, 
                                           ! the parameter describing the decrease rate 
                                           ! of methane production term 
                                           ! (due to "new" organics decomposition ~ NPP)
                                           ! with soil depth, m**(-1)

real(8), parameter :: r0_oldorg = 7.3d-7 !1.67d-7  ! the constant rate factor for "old" organics decomposition,  
                                                   ! according to formulation based on first-order kinetics, mol/(m**3*s)
real(8), parameter :: r0_oldorg_star = 4.5d-8 !1.67d-7  ! the constant rate factor for "old" organics decomposition,  
                                                   ! according to formulation based on first-order kinetics, mol/(kg*s)
real(8), parameter :: r0_oldorg2 = 7.3d-8          ! the constant rate factor for "old" organics decomposition, 
                                                   ! according to Michaelis-Menten based formulation, mol/(m**3*s)
real(8), parameter :: r0_oldorg2_star = 1.d-10 !2.d-10 ! the constant rate factor for "old" organics decomposition, 
                                                   ! according to Michaelis-Menten based formulation, mol/(kg*s)
real(8), parameter :: alpha_old_org = 2.d-1 !2.d-1 ! the constant of decomposition, year**(-1), value needs to be verified (!)
real(8), parameter :: k_oldorg = 3.d-1 ! the constant in the denominator of Michaelis-Menten equation when deriving 
                                       ! the equation for methane generation from old organics decomposition, kg/m**3
                                       ! varies from 0.1 to 0.3 (rough interval)
real(8), parameter :: V_oldorg = 0.2d-2 !0.2d-2 ! the constant in the numinator of Michaelis-Menten equation when deriving 
                                       ! the equation for methane generation from old organics decomposition, kg/m**3/year
                                       
real(8), parameter :: C0_oldorg = 18.  ! the density of organics sequestered under the talik, kg/m**3
real(8), parameter :: ice_trap_bubbl_meth_decr = 0.625 ! the fraction of methane that remains in bubbles
                                                       ! when they have been trapped by the ice cover;
                                                       ! the decrease of methane happens due to oxidation and 
                                                       ! the exchange with gases dissolved water
       
!*    =============================================================
!*               parameters for oxidation
!*             =============================

real(8), parameter :: vmax = 45./3600./1000. ! Michaelis-Menten, mol/(m**3*s) !(muM/s)
real(8), parameter :: rkm = 5./1000.         ! coefficients, mol/(m**3*s) !(muM)
real(8), parameter :: q10 = 2.               ! the observed value for oxidation, n/d

real(8), parameter :: vq_max = 6.d-4 ! maximum methane oxidation potential (10 deg C), 
                                     ! after Arah&Stephen (1998), Watson et al. (1997)
real(8), parameter :: temp0 = 283. ! reference temperature in Arrhenius equation, K
real(8), parameter :: k_ch4 = 0.44 ! Michaelis constant for methane in methane oxidation,
                                   ! after Arah&Stephen (1998), Nedwell&Watson (1995) 
real(8), parameter :: k_o2 = 0.33  ! Michaelis constant for oxygen in methane oxidation, 
                                   ! after Arah&Stephen (1998), Nedwell&Watson (1995)
real(8), parameter :: delta_Eq = 5.d+4 ! activation energy for methane oxidation, J/mol, 
                                       ! after Arah&Stephen (1998), Dunfield et al. (1993),
                                       ! Nedwell&Watson (1995)
       
!*	 =============================================================
!*                parameters for ebullition
!*              =============================

real(8), parameter :: rke = 1./3600.      ! a rate constant of the unit, 1/s
real(8), parameter :: cmin = 500./1000.   ! the concentration at which bubble formation occurs, mol/(m**3*s) !(muM)
real(8), parameter :: punveg = 0.         ! the percentage of unvegetated, bare soil
real(8), parameter :: cthresh = cmin*(1.+punveg/100.) !the threshold concentration for bubble formation
real(8), parameter :: rel_conc_ebul_crit = 0.4 !0.4, 0.45 the threshold relative concentration of methane
                                          ! for ebullition; 
                                          ! relative concentration = concentration divided by 
                                          ! maximal concentration according to Henry law at a given
                                          ! pressure
real(8), parameter :: meth_ebul_ice_intercept = 0.9 ! The portion of gas (methane) bubbles
                                                    ! intercepted by the ice cover during
                                                    ! a wintertime, n/d

!*	 =============================================================
!*            parameters for plant-mediated transport
!*          ===========================================

real(8), parameter :: rkp = 0.01/3600.    ! a rate constant of the unit 0.01/s
real(8), parameter :: tveg = 15.          ! a factor describing the quality of plant-mediated transport at a site, n/d
real(8), parameter :: pox = 0.5           ! a certain fraction of methane oxidized in
                                          ! the oxic zone around the roots of plants, n/d
real(8), parameter :: rlmin = 0.          ! the parameter for fgrow(t), n/d
real(8), parameter :: rl = 4.             ! the parameter for fgrow(t), n/d
real(8), parameter :: rlmax = 4.          ! rlmax=rlmin+rl, n/d


END MODULE METH_OXYG_CONSTANTS
